#include <fstream>
#include <iostream>
#include <iomanip>
using namespace std;

#include "plant_2DSOS_RDM.h"
// --------------------------------------------------
plant_2DSOS_RDM::plant_2DSOS_RDM(double idt,int iL,int iMode, double iE_D,double iE_S,double iE_N)
{
  // Options
  LatticeSize = iL;
  mode = iMode;
  assert(2*iMode<=iL);
  E_D = iE_D;
  E_S = iE_S;
  E_N = iE_N;

  // 2 Inputs: #0 = T, #1 = W
  int input_num = 2;
  int input_dim[2] = {1,1};
  // 4 outputs:
  //  Output #0~3: surface profile
  //  Output #4: h
  //  Output #5: r2
  //  Output #6: m2
  int output_num = 7;
  int output_dim[7];
  output_dim[0] = (mode+1)*(mode+1);
  output_dim[1] = (mode+1)*(mode+1);
  output_dim[2] = (mode+1)*(mode+1);
  output_dim[3] = (mode+1)*(mode+1);
  output_dim[4] = 1;
  output_dim[5] = 1;
  output_dim[6] = 1;
  block::initialize(input_num,input_dim,output_num,output_dim,idt);

  cov_z1 = outputs[0];
  cov_z2 = outputs[1];
  cov_z3 = outputs[2];
  cov_z4 = outputs[3];

  // States
  lattice = new int*[LatticeSize];
  lattice_stat =  new int*[LatticeSize];
  for(int i=0;i<LatticeSize;i++)
  {
    lattice[i] = new int[LatticeSize];
    lattice_stat[i] = new int[LatticeSize];
  }
  p_migr = new double[LatticeSize];
  p_desorp = new double[LatticeSize];

  ef_cos_table = new double*[mode+1];
  ef_sin_table = new double*[mode+1];
  int m,i;
  double SQRT_TWO = sqrt(2.0);
  double SQRT_PI = sqrt(pi);

  for(m=0;m<=mode;m++){
    ef_cos_table[m] = new double[LatticeSize];
    ef_sin_table[m] = new double[LatticeSize];
    for(i=0;i<LatticeSize;i++){
      if(m == 0){	
	ef_cos_table[m][i] = 1.0/SQRT_PI;
      }else{ 
	ef_cos_table[m][i] = SQRT_TWO/SQRT_PI*cos(2*m*pi/LatticeSize*i);
      }
      ef_sin_table[m][i] = SQRT_TWO/SQRT_PI*sin(2*m*pi/LatticeSize*i);
    }
  }
  reset();
}
// --------------------------------------------------
plant_2DSOS_RDM::~plant_2DSOS_RDM()
{
  int i;
  for(i=0;i<LatticeSize;i++)
  {
    delete []lattice[i];
    delete []lattice_stat[i];
  }
  delete []lattice;
  delete []lattice_stat;
  delete []p_migr;
  delete []p_desorp;
}
// --------------------------------------------------
void plant_2DSOS_RDM::reset()
{
  block::time = 0.0;

  // Inputs
  T = 0;
  W = 0;

  // States
  int x,y;
  for(y=0;y<LatticeSize;y++){
    for(x=0;x<LatticeSize;x++){
      lattice[y][x] = 0;
    }
  }
  for(y=0;y<LatticeSize;y++){
    for(x=0;x<LatticeSize;x++){
      init_local_stat(x,y);
    }
  }

  // Outputs
  int m,n,index;
  for(m=0;m<=mode;m++){
    for(n=0;n<=mode;n++){
      index = index2D_mode(m,n);
      cov_z1[index] = 0.0;
      cov_z2[index] = 0.0;
      cov_z3[index] = 0.0;
      cov_z4[index] = 0.0;
    }
  }
  h  = 0.0;
  r2 = 0.0;
  m2 = 0.0;
}
// ---------------------------------------------------
void plant_2DSOS_RDM::init_local_stat(int x,int y)
{
  lattice_stat[y][x]=(lattice[y][(x+LatticeSize-1)%LatticeSize]>=lattice[y][x])
    +(lattice[y][(x+LatticeSize+1)%LatticeSize]>=lattice[y][x])
    +(lattice[(y+LatticeSize-1)%LatticeSize][x]>=lattice[y][x])
    +(lattice[(y+LatticeSize+1)%LatticeSize][x]>=lattice[y][x]);
  return;
}
// ---------------------------------------------------
void plant_2DSOS_RDM::init_migr_rate()
{
  int i;
  double k_m;

  k_m=2*T*KB2H;

  for(i=0;i<4;i++)
  {
    MIGR_RT[i]=k_m*exp(-(E_S+i*E_N)*EV2KB/(T+0.0));
  }
  MIGR_RT[4]=0;
  return;
}
// ---------------------------------------------------
void plant_2DSOS_RDM::init_desorp_rate(){
  int i;
  double k_m;

  k_m=2*T*KB2H;

  for(i=0;i<5;i++)
  {
    //DESORP_RT[i]=k_m*exp(-(E_D+i*E_N)*EV2KB/(temperature+0.0));
    DESORP_RT[i]=0;
  }
  return;
}
// ---------------------------------------------------
void plant_2DSOS_RDM::init_p_migr()
{
  int x,y;

  migr_ttl=0.0;
  for(y=0;y<LatticeSize;y++)
  {
    p_migr[y]=0.0;
    for(x=0;x<LatticeSize;x++)
    {
      p_migr[y]+=MIGR_RT[lattice_stat[y][x]];

    }
    migr_ttl+=p_migr[y];
  }
  return;
}
// ---------------------------------------------------
void plant_2DSOS_RDM::init_p_desorp()
{
  int x,y;

  desorp_ttl=0.0;
  for(y=0;y<LatticeSize;y++)
  {
    p_desorp[y]=0.0;
    for(x=0;x<LatticeSize;x++)
    {
      p_desorp[y]+=DESORP_RT[lattice_stat[y][x]];

    }
    desorp_ttl+=p_desorp[y];
  }
  return;
}
// ---------------------------------------------------
void plant_2DSOS_RDM::update(double sysTime,double iT,double iW){
  if(time>sysTime)
    return;

  double rdm;
  while(time<=sysTime){
    if(iT!=T)
      set_temperature(iT);
    if(iW!=W)
      set_ads_rate(iW);

    rate_ttl = ads_ttl+migr_ttl+desorp_ttl;
    rdm = RandRand()*rate_ttl;
    if(rdm<ads_ttl)
      adsorption();
    else if(rdm<=migr_ttl+ads_ttl)
      migration();
    else
      desorption();

    time += OneMCTime(rate_ttl);
  }
  h2mode();
  outputs[4][0] = getH();
  outputs[5][0] = getR2();
  outputs[6][0] = getM2();
  
  return;
}
// --------------------------------------------------
void plant_2DSOS_RDM::update(double sysTime){
  update(sysTime,inputs[0][0],inputs[1][0]);
}
// --------------------------------------------------
void plant_2DSOS_RDM::getSnapshot(const string FileName)
{
  ofstream fout(FileName.c_str());
  fout << "LatticeSize  = " << LatticeSize << "\n";
  fout << "LatticeSize  = " << LatticeSize << "\n";
  fout << "Time = " << time << "\n";
  int x,y;
  for(y=0;y<LatticeSize;y++)
  {
    for(x=0;x<LatticeSize;x++)
    {
      fout << setw(5) << lattice[y][x];
    }
    fout << "\n";
  }
  fout.close();
}
// --------------------------------------------------
double plant_2DSOS_RDM::getM2()
{
  int x,y;
  double sum = 0;
  double m2;
  for(y=0;y<LatticeSize;y++)
  {
    for(x=0;x<LatticeSize-1;x++)
    {
      sum+=(lattice[y][x]-lattice[y][x+1])*(lattice[y][x]-lattice[y][x+1]);
    }
    sum+=(lattice[y][LatticeSize-1]-lattice[y][0])*(lattice[y][LatticeSize-1]-lattice[y][0]);
  }
  m2 = sum/LatticeSize/LatticeSize;
  return m2;
}
// --------------------------------------------------
double plant_2DSOS_RDM::getH()
{
  int x,y;
  double h_avg,sum = 0;
  for(y=0;y<LatticeSize;y++)
  {
    for(x=0;x<LatticeSize;x++)
    {
      sum+=lattice[y][x];
    }
  }
  h_avg = sum/LatticeSize/LatticeSize;
  return h_avg;
}
// --------------------------------------------------
double plant_2DSOS_RDM::getR2()
{
  int x,y;
  double R2;
  double sum=0;
  double h_avg;

  h_avg = getH();
  for(y=0;y<LatticeSize;y++)
  {
    for(x=0;x<LatticeSize;x++)
    {
      sum+=(lattice[y][x]-h_avg)*(lattice[y][x]-h_avg);
    }
  }
  R2 = sum/LatticeSize/LatticeSize;
  return R2;
}
// --------------------------------------------------
double plant_2DSOS_RDM::getTime()
{
  return time;
}
// --------------------------------------------------
void plant_2DSOS_RDM::set_temperature(double T_set)
{
  T = T_set;
  init_migr_rate();
  init_desorp_rate();
  init_p_migr();
  init_p_desorp();
  k_m2d = exp(-(E_D-E_S)*EV2KB/(T+0.0));
  return;
}
// --------------------------------------------------
void plant_2DSOS_RDM::set_ads_rate(double W_set)
{
  W=W_set;
  ads_ttl=W*LatticeSize*LatticeSize;
  return;
}
// --------------------------------------------------
void plant_2DSOS_RDM::add_remove_atom(int x, int y, int num)
{
  int xpos[5],ypos[5];
  double site_rate_mig, site_rate_desorp;
  int i;

  GetNeighbors(x,y,xpos,ypos);
  for(i=0;i<5;i++)
  {
    site_rate_mig = MIGR_RT[lattice_stat[ypos[i]][xpos[i]]];
    site_rate_desorp = DESORP_RT[lattice_stat[ypos[i]][xpos[i]]];
    migr_ttl  -= site_rate_mig;
    desorp_ttl -= site_rate_desorp;
    p_migr[ypos[i]] -= site_rate_mig;
    p_desorp[ypos[i]] -= site_rate_desorp;
  }

  if((lattice[y][x]==0) && (num==-1))
  {}
  else
  {
    lattice[y][x]+=num;
  }

  for(i=0;i<5;i++) 
  {
    init_local_stat(xpos[i],ypos[i]);
  }

  for(i=0;i<5;i++)
  {
    site_rate_mig=MIGR_RT[lattice_stat[ypos[i]][xpos[i]]];
    site_rate_desorp=DESORP_RT[lattice_stat[ypos[i]][xpos[i]]];
    migr_ttl+=site_rate_mig;
    desorp_ttl+=site_rate_desorp;
    p_migr[ypos[i]]+=site_rate_mig;
    p_desorp[ypos[i]]+=site_rate_desorp;
  }
  return;
}
// --------------------------------------------------	
void plant_2DSOS_RDM::GetNeighbors(int x, int y, int *xpos, int *ypos)
{
  xpos[0]=x;
  xpos[1]=x-1;
  xpos[2]=x;
  xpos[3]=x+1;
  xpos[4]=x;
  if(xpos[1]<0) xpos[1]+=LatticeSize;
  if(xpos[3]>=LatticeSize) xpos[3]-=LatticeSize;

  ypos[0]=y;
  ypos[1]=y;
  ypos[2]=y-1;
  ypos[3]=y;
  ypos[4]=y+1;
  if(ypos[2]<0) ypos[2]+=LatticeSize;
  if(ypos[4]>=LatticeSize) ypos[4]-=LatticeSize;

  return;
}
// --------------------------------------------------	
void plant_2DSOS_RDM::adsorption()
{
  add_remove_atom((int)(RandN1()*LatticeSize),(int)(RandN1()*LatticeSize),1);

  return;
}
// --------------------------------------------------	
void plant_2DSOS_RDM::migration()
{
  int x,y;
  int xpos[5],ypos[5];
  int i,n;
  int st[4][2];

  migr_site(&x,&y);
  if(lattice_stat[y][x]==4)
    return;

  GetNeighbors(x,y,xpos,ypos);

  // Count the number of neighbors with lower height
  for(i=1,n=0;i<5;i++)
  {
    if((lattice[y][x]-lattice[ypos[i]][xpos[i]])>0)
    {
      st[n][0]=ypos[i];
      st[n][1]=xpos[i];
      n++;
    }
  }

  if(n==0) return;

  n = (int)(RandN1()*n);
  add_remove_atom(x,y,-1);
  add_remove_atom(st[n][1],st[n][0],+1);

  return;
}
// --------------------------------------------------
void plant_2DSOS_RDM::desorption()
{
  int x,y;
  migr_site(&x,&y);
  add_remove_atom(x,y,-1);
  return;
}
// --------------------------------------------------
void plant_2DSOS_RDM::migr_site(int *migr_x, int *migr_y)
{
  int y,x;
  double rdm;

  rdm = RandN1()*migr_ttl;
  for(y=0;rdm>=0;y++)
    rdm-=p_migr[y];
  y--;
  rdm+=p_migr[y];
  for(x=0;rdm>=0;x++)
    rdm -= MIGR_RT[lattice_stat[y][x]];
  x--;

  *migr_x = x;
  *migr_y = y;

  /* FOR DEBUG
  int event[5];
  for(i=0;i<5;i++)
  {
  event[i] = 0;
  }
  if(lattice_stat[y][x]==0)
  event[0]++;
  else if(lattice_stat[y][x] == 1)
  event[1]++;
  else if(lattice_stat[y][x] == 2)
  event[2]++;
  else if(lattice_stat[y][x] == 3)
  event[3]++;
  else if(lattice_stat[y][x] == 4)
  event[4]++;
  */
  return;
}
// --------------------------------------------------
void plant_2DSOS_RDM::desorp_site(int *desorp_x, int *desorp_y)
{
  int y,x;
  double rdm;

  rdm = RandN1()*migr_ttl;
  for(y=0;rdm>=0;y++)
    rdm-=p_desorp[y];
  y--;
  rdm+=p_migr[y];
  for(x=0;rdm>=0;x++)
    rdm -= DESORP_RT[lattice_stat[y][x]];
  x--;

  *desorp_x = x;
  *desorp_y = y;

  return;
}
// --------------------------------------------------
double plant_2DSOS_RDM::OneMCTime(double r_d)
{
  double rd;

  rd=RandRand();
  while ((rd <= 0) || (rd >= 1))
    rd=RandRand();

  return -log(rd)/r_d;
}
// --------------------------------------------------
void plant_2DSOS_RDM::h2mode(){
  int m,n,i,j;
  int index;
  time_t time_begin;
  time_t time_end;
  double dx_square = (pi/LatticeSize)*(pi/LatticeSize);
  std::time(&time_begin);
  for(m=0;m<=mode;m++){
    for(n=0;n<=mode;n++){
      index = index2D_mode(m,n);
      cov_z1[index] = 0.0;
      cov_z2[index] = 0.0;
      cov_z3[index] = 0.0;
      cov_z4[index] = 0.0;
      for(i=0;i<LatticeSize;i++){
	for(j=0;j<LatticeSize;j++){
	  cov_z1[index]+=lattice[i][j]*phi1mn(m,n,i,j);
	  cov_z2[index]+=lattice[i][j]*phi2mn(m,n,i,j);
	  cov_z3[index]+=lattice[i][j]*phi3mn(m,n,i,j);
	  cov_z4[index]+=lattice[i][j]*phi4mn(m,n,i,j);
	}
      }
      cov_z1[index]*= dx_square;
      cov_z1[index] = cov_z1[index]*cov_z1[index];
      cov_z2[index]*= dx_square;
      cov_z2[index] = cov_z2[index]*cov_z2[index];
      cov_z3[index]*= dx_square;
      cov_z3[index] = cov_z3[index]*cov_z3[index];
      cov_z4[index]*= dx_square;
      cov_z4[index] = cov_z4[index]*cov_z4[index];
    }
  }
  std::time(&time_end);
  printf("h2mode() takes %f\n",std::difftime(time_end,time_begin));
}
/************************************************************************/
void plant_2DSOS_RDM::h2mode_fftw(){
  time_t time_begin;
  time_t time_end;
  std::time(&time_begin);

  std::time(&time_end);
  printf("h2mode() takes %f\n",std::difftime(time_end,time_begin));
}
///************************************************************************/
///* Choose and return an initial random seed based on the current time.
//Based on code by Lawrence Kirby <fred@genesis.demon.co.uk>.
//Usage: srand (time_seed ()); */
//unsigned time_seed (void)
//{
//  time_t timeval;	/* Current time. */
//  unsigned char *ptr;	/* Type punned pointed into timeval. */
//  unsigned seed;	/* Generated seed. */
//  size_t i;
//
//  timeval = time (NULL);
//  ptr = (unsigned char *) &timeval;
//
//  seed = 0;
//  for (i = 0; i < sizeof timeval; i++)
//    seed = seed * (0xff + 2u) + ptr[i];
//
//  return seed;
//}
///*****************************************************************************/
//double ranO(int *idum)
///*Returns a uniform random deviate between 0.0 and 1.0 using rand().*/
///*Set idum to any negative value to initialize or reinitialize the sequence*/
//{
//  static double y, maxran, v[131];
//  double dum;
//  static int iff=0;
//  int j;
//
//  if (*idum < 0 || iff == 0)
//  {
//    iff = 1;
//    maxran = RAND_MAX+0.0;
//    srand((unsigned)time(NULL));
//    *idum = 1;
//    for (j=0; j<=130; j++)
//    {
//      dum=rand()+0.0;
//    }
//    for (j=0; j<=130; j++)
//    {
//      v[j]=rand()+0.0;
//    }
//    y=rand()+0.0;
//  }
//  j= (int)(130.0*y/maxran);
//  if (j > 130 || j < 0)
//  {
//    printf("RanO: This cannot happen.\n");
//  }
//  y=v[j];
//  v[j]=rand();
//  return y/(maxran+1.0);
//}
//
///************************************************************************/
//double RandRand()
//{
//  double r_tot, r;
//  int i, n;
//  static int seed;
//
//  n=3;
//  r_tot=0.0;
//
//  for (i=0; i<n; i++)
//  {
//    r=ranO(&seed);
//    if (i > 0)
//    {
//      r/=pow((RAND_MAX+0.0), (i+0.0));
//    }
//    r_tot+=r;
//  }
//
//  return r_tot;
//}
//double RandN1()
//{
//  double rdm;
//
//  while((rdm=RandRand())>=1.0) {}
//
//  return rdm;
//}
/************************************************************************/